%%
close all
clear;
N = 2;
L = 4;
length1 =L*N;
bitlist =randi([0 1], 1,length1);
numslist = bi2de(reshape(bitlist, length(bitlist)/N, N));
% TS_A=qammod(numslist, 2^N, 'UnitAveragePower', true).';
TS_A = pskmod(numslist, 2^N, pi/4).';
TS_A = repmat(TS_A,1,8);
fs = 5000;
Nfft = 64;

%f移到中点
f = (-Nfft/2:Nfft/2-1)*fs/Nfft;
% fs/Nfft*8;
figure;
subplot(2,1,1);
plot(f,abs(fftshift(fft(TS_A,Nfft))));
title("频谱密度")
%给TS_A增加频偏
TS_A_offset = TS_A.*exp(-1i*2*pi*0.003*(1:L*8));
subplot(2,1,2);
plot(f,abs(fftshift(fft(TS_A_offset,Nfft))),'r');,
title("增加频偏后的频谱密度")
scatterplot(TS_A_offset);

[c,lags] = xcorr(abs(fftshift(fft(TS_A,Nfft))),abs(fftshift(fft(TS_A_offset,Nfft))));
figure;stem(lags,c);title("互相关函数");
[~,flag] = max(c);
offset = lags(flag)/Nfft;
fprintf('粗频偏估计是: %d\n', offset);
% end

%%
TS_A_NEW = zeros(1,length(TS_A)*2);
TS_A_NEW(1:2:end)=TS_A;


beta=0.125; %%滚降系数
NumSym=length(TS_A_NEW);
HRRC=zeros(1,NumSym);
HRRC([1:floor((1-beta)*NumSym/2)+1,end-floor((1-beta)*NumSym/2+1):end])=1;
index=floor((1-beta)*NumSym/2)+1:floor((1+beta)*NumSym/2)-1;
HRRC(index+1)=(1+cos(pi*2/beta*((index-floor((1-beta)*NumSym/2))/NumSym/2)))/2;
HRRC=sqrt(HRRC);
HRRC(end-index+1)=HRRC(index+1);
TS_A_NEW =ifft(fft(TS_A_NEW).*HRRC);


nodata=TS_A_NEW;
rand_nums = 0.6;%%引入时钟误差
[dataRTM]=sqretiming_noise(nodata,rand_nums);    
Pu_out_noise=dataRTM/mean((abs(dataRTM)));
scatterplot(Pu_out_noise);
title("引入时钟误差后的星座图")
[dataRTM,restore_val] = Sqretiming1(Pu_out_noise);
scatterplot(dataRTM);


Xk = fftshift(fft(Pu_out_noise,Nfft));
Xk_conj = conj(circshift(Xk,-Nfft/2+1));
total = Xk.*Xk_conj;
total2 = 0;
for ii = 1:Nfft/2
    total2=total2+total(ii);
end

tao = imag(total2);

%%
clear;
N = 2;
length1 =10*N;
bitlist =randi([0 1], 1,length1);
numslist = bi2de(reshape(bitlist, length(bitlist)/N, N));
TS_B=qammod(numslist, 2^N, 'UnitAveragePower', true)';
%TS_B添加伪随机噪声
TS_B_noise = awgn(TS_B,20);
TS_B_total = [TS_B,TS_B_noise];
TS_B_total = repmat(TS_B_total,1,3);

[c,lags] = xcorr(TS_B_total,TS_B);
figure;stem(lags,c);title("互相关函数");
[~,flag] = max(c);
offset = lags(flag);
disp("细频偏",offset)

%%
TS_B_offset = TS_B_total.*exp(-1i*2*pi*0.03*(1:length1*3));
var1 = angle(TS_B_offset(11:20)*TS_B_offset(1:10)');
offset = var1/10/(2*pi);
scatterplot(TS_B_offset);
title("未补偿细频偏")
res = TS_B_offset.*exp(-1i*2*pi*(offset)*(1:length1*3));
scatterplot(res);
title("补偿细频偏")

